#include <iostream>
#include <math.h>

using std::cout;
using std::endl;


extern "C" { // Derek's smith code

  void def_universe_(const float& om_m, const float& om_v,
      const float& om_b, const float& h, const float& sig_8,
      const float& w, const float& wa);
  float smithp3d_(const float& k, const float& z);
}

int main()
{
  const float om_m = 0.27;
  const float om_v = 0.73;
  const float om_b = 0.046;
  const float h = 0.72;
  const float sig_8 = 0.9;
  const float w = -1.;
  const float wa = 0.;
  const float w2 = -0.8;
  const float wa2 = -0.6;
  const float z = 1.;

  const float kfactor = exp(log(1.e5)/99);

  float k[100], light[100], dark[100];

  k[0] = 0.01;
  for(int i=1;i<100;i++) k[i] = k[i-1]*kfactor;

  def_universe_(om_m,om_v,om_b,h,sig_8,w,wa);
  for(int i=0;i<100;i++) light[i] = smithp3d_(k[i],z);
  def_universe_(om_m,om_v,om_b,h,sig_8,w2,wa2);
  for(int i=0;i<100;i++) dark[i] = smithp3d_(k[i],z);

  cout << "# k(h/Mpc) P_Smith P_Smith_DE\n";
  for(int i=0;i<100;i++) 
    cout << k[i] << "  " << light[i] << "  " << dark[i] << endl;

  return 0;
}
    
